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The performance of a state-of-the-art continuum damage mechanics model for 
interlaminar damage, coupled with a cohesive zone model for delamination is examined for 
failure prediction of quasi-isotropic open-hole tension laminates. Limitations of continuum 
representations of intra-ply damage and the effect of mesh orientation on the analysis 
predictions are discussed. It is shown that accurate prediction of matrix crack paths and 
stress redistribution after cracking requires a mesh aligned with the fiber orientation. Based 
on these results, an aligned mesh is proposed for analysis of the open-hole tension specimens 
consisting of different meshes within the individual plies, such that the element edges are 
aligned with the ply fiber direction. The modeling approach is assessed by comparison of 
analysis predictions to experimental data for specimen configurations in which failure is 
dominated by complex interactions between matrix cracks and delaminations. It is shown 
that the different failure mechanisms observed in the tests are well predicted. In addition, 
the modeling approach is demonstrated to predict proper trends in the effect of scaling on 
strength and failure mechanisms of quasi-isotropic open-hole tension laminates. 


I. Introduction 

A luminum structures for transport category aircraft have been in service for many years, and the design 
processes for these structures are well established. Design and certification procedures have been developed to 
ensure safe operations in the event of local damage or cracks in the structure. These design and certification 
procedures are based on the damage tolerance design philosophy that requires that any propagating crack in a 
structure be arrested by the next contiguous primary structural element. This design philosophy is based on linear 
structural analysis methods to determine the stresses in the aircraft structure and on linear elastic fracture mechanics 
to determine the crack growth characteristics. Extensive structural tests are usually relied upon to validate this linear 
structural design approach. 

The introduction of laminated composite materials to aircraft structural applications has led to the development 
of similar linear elastic, fracture mechanics-based analysis and damage tolerant design approaches for composite 
structures. Numerous analytical techniques have been developed for predicting the notched strength of flat 
composite laminates subjected to uniaxial loads (e.g. Refs. 1-7). A comprehensive review of these models has been 
provided by Awerbuch and Madhukar. 7 A large majority of the techniques that have been developed are semi- 
empirical, strength-based or fracture mechanics-based theories. Their development was motivated by the need to 
provide simple and efficient methods for predicting the notched strength of composite structures. Most of these 
models have demonstrated success in predicting laminate strength in the presence of holes or notches, provided that 
the parameters in the model are properly determined. In general, the parameters in the models are assumed to be 
independent of specimen geometry. However, the parameters are strongly dependent on the material system and 
laminate configuration and, consequently, must be determined for each new material system and laminate stacking 
sequence. In addition, none of the semi-empirical fracture models address the characteristics of the crack- tip damage 
zone, which are laminate and material dependent, or the complex interaction of micro- and macro-failures associated 
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with the crack-extension process. Instead, the crack-tip damage zone is simulated as some “effective” notch-tip 
damage zone that is assumed to grow in a self-similar manner. In many cases, self-similar crack growth typical of a 
similar metal component is not observed. 

The potential benefit of augmenting current design and certification processes with high-fidelity analysis 
methods to reduce empiricism and extensive structural testing is well recognized. Quantification of the structural 
fracture resistance under various loading conditions is fundamental for evaluating the residual strength and life of 
aircraft structures. This requirement has led to the development of progressive damage models to predict damage 
initiation, damage growth and ultimate strength of composite laminates. Significant progress has been achieved in 
developing numerical approaches and underlying constitutive models for initiation and propagation of specific 
damage modes. Intra-ply damage modes have been investigated primarily within the framework of Continuum 
Damage Mechanics (CDM), 8 ’ 9 while delamination has been studied extensively using interface fracture modeling 
techniques such as cohesive zone models 10 and virtual crack closure techniques (VCCT). 10,11 Despite these advances 
in progressive damage modeling, recent studies 12 " 15 have suggested that CDM models for ply failure, coupled with 
cohesive zone models for delamination have difficulty in accurately representing laminate failure sequences that 
exhibit strong coupling between transverse matrix cracking and delamination. In the CDM methodology a crack or 
displacement discontinuity is replaced with local volumetric stiffness degradation. As the crack opens, the stiffness 
degradation is such that stress should not be transferred across the crack faces. Iarve 15 showed that this is not 
necessarily the case, and because of the spurious stress transfer, local stress redistribution is not accurately predicted 
in the vicinity of damage. The stress locking or spurious stress transfer observed in these studies, is mesh orientation 
dependent and may be alleviated by aligning the element edges with the crack direction. 16 ' 19 

The present paper presents a study of the performance of a continuum damage mechanics model for intralaminar 
damage, coupled with a cohesive zone model for delamination, for failure prediction of quasi-isotropic open-hole 
tension specimens. This specimen configuration has been studied extensively, both experimentally and 
computationally, 20 " 24 and represents a significant challenge for progressive damage analysis for laminate failure that 
is dominated by complex interactions between transverse matrix cracks and delaminations. First, limitations of state- 
of-the-art continuum representations of intra-ply damage, and mesh-objectivity are discussed. Particular emphasis is 
given to the effect of mesh orientation on the analysis predictions and the effect of mesh orientation on stress 
locking. Analysis predictions of a unidirectional open-hole tension specimen obtained with two meshes, a radial 
mesh, and a mesh with element edges aligned with the fiber orientation are presented to evaluate stress locking. 
Results of these analyses are then used as motivation for the development of a finite element model for the quasi- 
isotropic laminates where individual plies are modeled with meshes that have element edges aligned with the fiber 
direction. Analysis predictions obtained with this mesh are then compared with experimental data to assess the 
modeling approach. 


II. Limitations of CDM Models 

To predict the ultimate strength of composite laminates and structures, an accurate numerical representation of 
all damage modes and their interactions is required. Some of the most complex damage models available 24 " 27 rely on 
CDM models to represent the intralaminar damage modes (e.g., transverse matrix cracking and fiber failure) and use 
cohesive zone models to capture delamination between ply interfaces. However, despite advances in progressive 
damage modeling, recent studies (e.g., van der Meer [13]) indicate that CDM models coupled with cohesive zone 
models may not always represent laminate failure sequences properly. These deficiencies are particularly evident 
when the observed fracture mode exhibits matrix splitting and pullouts 20 or when the fracture is characterized by a 
strong coupling between transverse matrix cracking and delamination. 13 

The deficiencies of the predictive capabilities can be attributed to several factors, including the incorrect 
prediction of the damage zone size normal to the fracture direction when using crack-band models, and the inability 
of local CDM models to reliably predict matrix crack paths. These limitations are mostly due to the fact that CDM 
models are usually implemented as “local” rather than “non-local” models 28 i.e., the evolution of damage in a local 
CDM model is evaluated at individual integration points without consideration of the state of damage at neighboring 
locations, and the inability of standard finite elements to represent localized shear bands. The following discussion 
pertains to local implementations of CDM models, since non-local damage models are less widely used due to the 
difficulty in implementing them within the finite element method. 

The premise of the crack-band approach for regularizing CDM models is that damage localizes into a band with 
a width equivalent to the element dimension. If the element size is smaller than the damage process zone, the crack- 
band approach may not predict correctly the width of the damage zone nor the local stress field. Consequently, the 


2 

American Institute of Aeronautics and Astronautics 



stress redistribution resulting from damage development may be inaccurately predicted and can potentially result in 
inaccurate representation of damage mode interactions and failure sequences. 

As a result of homogenization and damage localization, CDM models have difficulties predicting crack paths. 
Since homogenization eliminates the distinction between fibers and matrix, a CDM model cannot distinguish 
between cracks that propagate along fiber directions from those that cross fibers. 13 In CDM models implemented 
with damage localization, the damage state at any integration point in the model depends only on the stress field at 
that point rather than the damage state of neighboring points. Therefore, the direction of damage evolution is driven 
only by the instantaneous local stress distribution. In other words, the local direction of cracking may be predicted 
correctly by the failure criteria, but, if the morphology of the material is not properly accounted for in the damage 
model, the sequence of failures that eventually defines the path of a crack at a macroscopic level may be predicted 
incorrectly. 
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(a) Propagation of shear damage along fiber (b) Propagation across fiber direction 


Figure 1. Idealized propagation of shear damage (adapted from [13]). 

The potential inability of CDM models to determine the correct direction of propagation is particularly evident 
when the stress field is dominated by shear. Consider two different plies in a laminate with a notch and subjected to 
shear, as shown in Fig. 1. In both situations, the stress level required to initiate matrix micro cracks is correctly 
predicted by the failure criterion. Furthermore, both situations would result in an identical sequence of failures, since 
the stress field is identical. However, it is clearly not the same to propagate a crack in a sequence of linked micro 
cracks (Fig. la) as it is to propagate a crack band across fibers (Fig. lb). Matrix cracking in a shear band running 
parallel to the fibers is a relatively brittle failure mechanism, whereas matrix cracking normal to the fibers produces 
a damage band that requires much more work to propagate. 

The sensitivity of CDM predictions to the finite element mesh orientation also contributes to the difficulty in 
predicting the crack path. Although the objectivity of the solution with respect to element size is addressed with the 
crack-band approach 29 the predicted damage may be dependent on mesh orientation and element shape. When 
strain-softening constitutive models are used in a finite element simulation, damage tends to propagate along 
preferred directions, consisting of either element edges or element diagonals. A demonstration of the sensitivity of 
simulation results to mesh orientation is provided in Fig. 2 for a unidirectional compact tension specimen with fibers 
oriented at 90° to the load direction. Results are presented for simulations obtained with a mesh oriented parallel to 
the fiber direction (Fig. 2a) and with an inclined mesh in front of the crack tip (Fig. 2b). The crack should propagate 
along the fiber direction. However, the results show directional bias, and the simulated crack band propagates in the 
direction of the element alignment. 

The tendency for damage to localize along mesh lines can be partially attributed to shear locking. 17 In the CDM 
methodology, a crack or displacement discontinuity is represented by a degradation of the corresponding terms in 
the constitutive stiffness. As the crack opens, the stiffness degradation is such that stress should not be transferred 
across the crack faces. However, unloading may not occur due to in-plane shear locking. Shear locking here refers to 
inappropriate shear stress transfer across a widely open smeared crack, which occurs when an element cannot shear 
without inducing tensile strains. In a ply with orthotropic properties, a matrix crack should be represented by setting 
the transverse shear modulus G 12 and the transverse Young’s modulus E 22 to zero. However, quadrilateral elements 
have been shown to exhibit coupling between the y 12 and e n strains unless the element edges are aligned with the 
softening band or are oriented at 45 degrees to the band. 17 Furthermore, the tendency of the element to lock is 
dependent on the order of integration of the element: fully integrated elements are more susceptible to pathological 
in-plane shear locking than reduced-integration elements. 
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(a) Mesh aligned with crack direction 


(b) Mesh inclined to crack direction 


Figure 2. Effect of mesh orientation on crack path in a unidirectional compact tension specimen. 


Shear stress transfer across an open smeared crack caused by shear locking can result in inaccurate prediction of 
stress redistribution after damage development. Iarve 15 demonstrated shear locking for a simple case of a 
unidirectional [0 8 ] graphite-epoxy open-hole tension specimen. Experiments show that this specimen exhibits 
splitting cracks parallel to the load and tangent to the hole. Splitting near holes and notches reduces the stress 
concentration, so predicting their effect is essential to obtain the ultimate strength with any accuracy. Iarve ’s 
progressive damage analyses were conducted using a radial-type mesh pattern with different levels of mesh 
refinement, where a radial mesh consists of the typical pattern of elements radiating from a circular hole and ending 
at a rectangular boundary. It was demonstrated that it is not possible to predict the stress relaxation with a radial 
mesh due to longitudinal splitting unless the fiber-direction modulus, E n , is also set to zero. 



(a) Radial Mesh (b) Aligned Mesh 


Figure 3. Splitting damage predicted using a radial mesh and an aligned mesh. 

When damage localizes and a fracture path is known a priori , a relatively simple approach to circumvent some 
of the limitations of CDM noted above consists of aligning the mesh with the direction of fracture. 19 Mesh alignment 
can be used to force a matrix crack to localize along the fiber direction. The benefit of mesh alignment is 
demonstrated below for an open-hole tension specimen. Predictions were obtained for an [Og] IM7-8552 laminate 
using the continuum damage mechanics model provided within the ABAQUS® finite element code. 32 Each ply was 
modeled using quadrilateral reduced integration shell elements, S4R, to minimize the tendency for locking. Analyses 
were obtained using a radial mesh and an aligned mesh, as shown in Fig. 3. The loading direction and fiber direction 
are parallel to the A-axis. 

The damage zones predicted in the region of the hole at the peak load, using the radial and aligned meshes are 
shown in Fig. 3. The results indicate that both models predict longitudinal splitting tangent to the hole. The load 
versus end-displacement response obtained with the two models is shown in Fig. 4a. The radial mesh model predicts 
failure of the specimen at 77 kips compared to the peak load of 116 kips predicted by the model with the aligned 
mesh. The peak load predicted with the radial mesh was significantly lower than the predicted peak load with the 
aligned mesh. This is probably due to the inability of rotated elements with the radial mesh to represent shearing 
along the axial split. The axial stress, an, normalized by the axial stress at the loading end of specimen, a x , obtained 
from the radial mesh and the aligned mesh at 77kips is shown as a function of distance from the hole edge (Y- 
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End displacement, in. 


(a) Load versus end-displacement 



Distance from hole, in. 

(b) Normalized axial stress distribution 


Figure 4. Effect of mesh type on predicted ultimate failure and fiber-stress relaxation due to splitting in a 
[0] 8 laminate. 


direction) in Fig. 4b to demonstrate effect of the mesh orientation on the stress distribution at the hole. The results 
obtained with the aligned mesh indicate that the stress concentration at the hole is nearly eliminated as a 
consequence of the split, as expected. The results obtained using the radial mesh show insufficient stress relaxation 
at the edge of the hole, indicating that damaged elements transfer shear load across the damaged band even though 
the matrix-dependent moduli G 12 and E 22 have been set to zero. 

III. Problem Description and Analysis Procedure 

A. Open-Hole Tension Test - Experimental Summary 

The open-hole tension (OHT) specimen was recently studied, both experimentally and computationally by 
Hallett and Green et al. 20,21 A schematic of the specimens considered is shown in Fig. 5. Specimens were fabricated 
from IM7/8552 graphite-epoxy and had a quasi-isotropic stacking sequence. Three-dimensional scaling, with ply- 
level and sublaminate-level approaches to through-the-thickness scaling, was examined to determine size effects in 
laminate failure. Ply-level scaled specimens are considered in the present study. Ply-level scaled specimens were 
fabricated by blocking multiple plies with the same orientation together, thereby increasing the effective ply 
thickness. Four thicknesses were considered with the stacking sequence [45 m /90 m /-45J0 m ] s and with m = [1, 2, 4, 8]. 
The plate width to hole diameter ratio, w/d was kept constant at a value of 5, and the length to hole diameter ratio, 
l/d, was equal to 20. 

The experimental data indicate that the failure progression and the mode of laminate failure are dependent on the 
laminate layup (stacking sequence and ply thicknesses) and that the scaling of strength with respect to size is 
dependent on the type of failure. Sub-critical damage that develops at relatively low loads prior to the ultimate 



I/d = 20 

Figure 5. Specimen configuration. 20 
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failure of the OHT specimens was determined to have a significant impact on the ultimate failure load and the 
failure mode. The progressive damage development in the OHT specimens was reported in detail by Hallet et al., 21 
and is summarized here. For all geometries, there were several damage processes that evolved in a coupled manner 
before final failure. The first failure observed in the tests was matrix cracking in the 90° plies at the free edge of the 
hole. These cracks were followed by the development of matrix cracks in the off-axis plies and in the 0° plies at the 
free edge of the hole, typically at tangent locations. Dominant cracks or splits developed in the surface 45° plies and 
in the 0° plies. The cracks in the -45° and 90° plies were more distributed and developed along the 45° and 0° splits. 
With continued loading, a triangular delamination developed at the hole edge between the surface 45° and the 
adjacent 90° ply. The splits also increased in length with continued loading until the split in the 45° surface plies 
reached the free edge. At that point, another triangular delamination region formed at the specimen free edge 
between the 45° split and a crack in the 90° ply. The damage at the specimen free edge then linked up with the 
damage at the hole edge, and propagated through the thickness of the laminate, until reaching the 0° ply. 

The extent of the sub-critical damage development described above is dependent on the laminate geometry and 
significantly influences the laminate failure mode. Three laminate failure modes were identified and were classified 
as either brittle, pull-out or delamination. For small values of blocked plies, the brittle failure mode was observed, in 
which little sub-critical damage develops and is dominated by fiber failure. For increasing numbers of blocked plies, 
a transition to a pull-out failure mode and then to a delamination failure mode was observed. The pull-out failure 
mode is characterized by sub-critical damage that propagates across the width of the specimen, allowing the off-axis 
plies to separate from each other, without fiber failure. The delamination failure mode is characterized by significant 
sub-critical damage development in the form of transverse matrix cracks in the individual plies and delaminations 
that link-up prior to ultimate failure. The three failure modes are shown in Fig. 6. 



(a) Brittle failure mode (b) Pull-out failure mode (c) Delamination failure mode 

Figure 6. Open-hole tension specimen failure modes. 20 


B. Progressive Damage Analysis 

Progressive damage analysis for the open hole tension specimens was conducted in ABAQUS explicit using the 
continuum damage mechanics model implemented in ABAQUS 30 to represent the intralaminar damage modes (e.g., 
transverse matrix cracking and fiber failure), and cohesive elements to capture delamination at ply interfaces. The 
intralaminar and interlaminar damage models implemented in ABAQUS are briefly described below. 

1. Intralaminar Damage 

The continuum damage mechanics model implemented in ABAQUS to simulate intralaminar damage is 
available for elements that have a plane stress formulation, and is intended for the prediction of elastic-brittle 
materials that show no appreciable plastic deformation before failure. Four failure modes including fiber tension, 
fiber compression, matrix tension, and matrix compression are considered and are represented separately. Initiation 
of damage, which refers to the onset of damage at a material point is based on Hashin’s theory. 31 The Hashin 
Criteria for the four different failure modes are described in Equations (1) - (4) as follows: 


Fiber tension (a n >0) 


F‘ 
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X) 
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Fiber compression (a x , < 0) 
Matrix tension (a 22 ^ 0) 

Matrix compression (a 22 < 0) 



( 2 ) 

( 3 ) 

( 4 ) 


In the above equations, a,y are the components of the effective stress tensor, and X T and X c are the longitudinal 
tensile and compressive strengths, Y r and Y° are the tensile and compressive strengths in the matrix direction, and S L 
and S T denote the longitudinal and transverse shear strengths. The coefficient a in Equation (1) determines the 
contribution of the shear stress to the initiation of fiber tensile failure. Once a damage initiation function is satisfied, 
the associated damage variable is different from zero and further loading will cause degradation of the material 
stiffness coefficients. The stiffness matrix of a damaged ply is defined as 

<J = C d S ( 5 ) 

(l ~df^)E n {\-d X\-d m) V 2\E 22 0 

^ d= ~^ (l - ^/ ) i^—d m ) V \2^22 ) ^22 0 (6) 

L 0 0 (wJg i2 d 

where D = d f )(l- d m )v l2 v 2l d f and df is the damage variable associated with fiber fracture, and d m , and d s , 

are damage variables associated with matrix failure. The damage evolution laws for the damage variables are 
defined in terms of the fracture energy dissipated during the damage process. Methods for determining the fracture 
energies are provided in Ref [32]. The damage evolution laws need to ensure that the computed energy dissipated is 
independent of the mesh. To do so, the energy dissipated for each damage mechanism is regularized using the crack- 
band model. 29 

2. Interlaminar Damage 

Interlaminar damage is simulated by placing ABAQUS cohesive elements, COH3D8, at potential delamination 
interfaces. The constitutive response of the cohesive elements in single-mode loading is defined by a bi-linear 
traction- separation law. Initially the crack is assumed to be elastic, and the crack closing forces are related to the 
interfacial displacement jump, 8 coh by a high penalty stiffness K. If the interfacial displacement jump exceeds a 
critical value, 8 h damage is assumed to have initiated, and the stiffness of the cohesive element is degraded. The 
crack closing forces are assumed to soften linearly such that the area under the traction-displacement curve is equal 
to the fracture toughness, G c . Complete separation is achieved when the displacement jump exceeds Sf. For single- 
mode loading, the bi-linear cohesive law can be expressed as: 

f KS, S<S l 

a = \(\-d*)KS, S t <S<S f (7) 

where the damage variable d* is a function of the displacement jump and accounts for the reduction in the load- 
carrying ability of the material as a result of damage. The damage variable d* has a value of zero when the interface 
is undamaged, and a value of one when the interface is fully fractured. 

For crack growth under mixed-mode loading, ABAQUS provides options for defining the point of damage 
initiation and the critical energy release rate for damage evolution. In all of the analysis presented herein, the onset 
of delamination is determined based on the quadratic interaction criterion 
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Table 1. IM7/8552 Material Properties 


E 11 

E 22 

G ,2 

V 12 

X, 

X c 

Y, 

Y c 

S 

(msi) 

(msi) 

(msi) 


(msi) 

(msi) 

(msi) 

(msi) 

(msi) 

24.85 

1.316 

0.765 

0.32 

0.3373 

0.1740 

0.0090 

0.02897 

0.0134 



where the Macaulay bracket < > indicates that a compressive normal stress <j zz does not contribute to the damage 
initiation, and T and S are the interfacial normal and shear strengths, respectively. The shear strengths in the two 
orthogonal directions are assumed to be equal. The critical fracture energy is determined from the Benzeggaagh and 
Kenane (BK) criterion 33 


G C = G ic+( G nc +G ic) 


J Shear 

v G t ) 


( 9 ) 


where 77 is the BK material parameter, G shear = G n + G IU , G T = Gj + G n + G ni , and G IC and G IIC are the critical 
fracture toughness values for pure Mode I and Mode II fracture, respectively. 

C. Finite Element Models 

Analyses were conducted for a subset of the ply-level scaled specimens tested by Winsom et al. and include 
specimens with d =1/8, 1/4, and 1/2 inches and with m = 1,2 and 4. Mechanical properties for the IM7/8552 
graphite-epoxy material are taken from Camahno , 24 and are provided in Table 1. Fracture toughness parameters for 
the fiber and matrix failure modes are provided in Table 2 . 24 The damage parameters provided for the matrix are 
used for both matrix damage and for interlaminar damage. To account for the in-situ effect on the intralaminar 
damage, the in-situ transverse tensile and shear strengths are determined using the unidirectional strength properties 
and the formulas proposed by Camanho et al . 34 The in-situ strengths take into account the ply thickness and the 
difference in constraint effects on an inner or outer ply. In-situ strength properties used in the analyses are provided 
in Table 3. In addition, the initial elastic stiffness for the interface cohesive law is taken to be K = 8.0 x 10 8 psi/in. 

Table 2. IM7/ 8552 Fracture Energies for Longitudinal and Transverse Frac ture (lb-in/in 2 ) 


Gic+ (fiber) 

Gic_ (fiber) 

G IC (matrix) 

G n c (matrix) 

461.5 

601.9 

1.571 

4.462 


Finite element models were developed to simulate coupled interlaminar and intralaminar damage progression. 
The thickness of each model was represented by seven sublaminates, with a cohesive element layer of zero thickness 
between them. Each sublaminate represents a group of blocked plies in a laminate. Two mesh types were used for 
the in-plane discretization to demonstrate the effect of mesh orientation on the predicted development of damage 


Table 3. Calculated In-situ Strengths for IM7-8552 



Orientation, thickness (in) 

Y t 

(psi) 

S 

(psi) 

Surface ply 

45°, 0.02 

7.31 x 10 3 

9.71 x 10 3 

Inner ply 

907-45°, 0.02 

1.43 x 10 4 

1.66 xlO 4 

Inner ply 

0°, 0.04 

1.43 xlO 4 

1.66 x 10 4 
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(a) In-plane discretization (b) Through-thickness discretization 

Figure 7. Finite element model, radial mesh. 

and specimen failure. The first mesh, referred to as a radial mesh, is shown in Fig. 7. The radial mesh was generated 
by aligning element edges along radial lines around the hole, and all plies were modeled with the same mesh. The 
cohesive elements between the plies were connected to the structural ply elements by coincident nodes. A second 
mesh, referred to as the aligned mesh, was also used. The aligned mesh consisted of different meshes within the 
individual plies, such that the element edges aligned with the ply fiber direction. This modeling approach is 
proposed on the basis that transverse matrix cracks in individual plies propagate along the fiber direction. The 
aligned mesh can potentially eliminate spurious stress transfer across crack faces that would otherwise occur, as has 
been demonstrated by Jirasek. 16-18 Individual finite element discretizations for each ply, and the cohesive element 
layer discretization are shown in Fig. 8 for the aligned mesh. Tie constraints were used to connect the ply structural 
elements with the cohesive layer elements. In-plane mesh refinement was used to provide elements with side lengths 
approximately equal to 0.01 inches (approximately 0.3 mm) in a region near the hole for the specimens with the 
smallest in-plane dimensions. This mesh density has been recommended to accurately represent the process zone for 
simulating delamination propagation 10 and for modeling progressive intralaminar damage. For the larger specimens, 
the elements were scaled by the in-plane scale factor, giving element edge lengths near the hole equal to 0.02 inches 
and 0.04 inches, for the specimens with !4-inch and ^-inch diameter holes, respectively. For these specimens the 
element edge lengths are larger than recommended, but this mesh density was retained to keep the computational 
requirements manageable. A detailed study of the effect of element size on the predictions has not been conducted, 
and will be addressed in future studies. 

Each sublaminate or blocked-ply group was modeled using the ABAQUS quadrilateral continuum shell element, 
SC8R. The SC8R continuum shell element is an 8-node, quadrilateral, first-order interpolation shell element with 
reduced integration and 3 degrees of freedom per node. Blocked plies were numerically represented by a single 
element through the thickness of the ply block. The cohesive layers between plies were modeled using the 
ABAQUS COH3D8 zero-thickness cohesive elements. The boundary conditions at the left end of the specimen were 
fully clamped ( u , v, w = 0). At the right end of the specimen, the transverse displacements, v and w were constrained 
and a uniform axial displacement, u, was prescribed. In the explicit analysis, a final displacement, and a simulation 



45° layer mesh 



-45° layer mesh 





Figure 8. Finite element model, aligned mesh. 
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Figure 9. Predicted load versus displacement response obtained with the radial mesh and the aligned 
mesh, d = 1/8-inch, m = 4. 


period are prescribed, to give a displacement loading rate. For computational efficiency, the displacement rate was 
between 2.0 (in./sec.) and 3.5 (in./sec.) in all analyses. A mass density of 1.5 x 10" 3 lb/in 3 was used for the ply 
structural elements, and the cohesive elements were assumed to be weightless. During the simulations, the kinetic 
energy was monitored to ensure that the kinetic energy did not exceed 5% of the total strain energy to ensure that 
significant dynamic effects were not introduced into the simulations by using relatively short total loading times. 

IV. Results 

Numerical results of the progressive damage analysis of the quasi-isotropic open-hole specimens are presented in 
this section. First, damage evolution and ultimate failure results predicted with the radial and aligned meshes are 
presented and compared with experimental observations. Second, the predicted effect of in-plane (2D) scaling and 
through-thickness scaling (ID) on the laminate strength and failure mode is compared with experimental data. 
Discrepancies between the analysis predictions and the experimental results are described, and possible explanations 
for the discrepancies are provided. 

A. Radial Mesh and Aligned Mesh Predictions, d = 1/8-inch, m = 4 

Predicted load versus deflection responses obtained with the radial mesh and the aligned mesh for the specimen 
with a 1/8-inch diameter hole and m = 4 are shown in Fig. 9. Test results obtained by Green, Wisnom and 
Hallett, 20,21 indicate that this specimen configuration fails in the delamination failure mode. For this mode, two 
failure stress levels were reported: an average failure stress and an ultimate failure stress. The average failure stress 
was defined as the stress level at which there was a significant load drop (5% or greater) in the load-displacement 
response. Table 4 shows the average failure and ultimate failure stress results obtained with the two meshes, along 
with the experimental data. As indicated in Fig. 9 and Table 4, results obtained with the aligned mesh are consistent 
with the experimental observations. A significant load drop is observed in the load versus displacement response 

Table 4. Experimental and Predicted Average and Ultimate Failure Stresses for Specimen with d = 1/8- 


inch, m = 4 


Experimental Results 

Radial Mesh 

Aligned Mesh 

Average Failure Stress, ksi 

39.9 

49.2 

48.2 

Ultimate Failure Stress, ksi 

65.3 

49.2 

60.7 
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Figure 10. Matrix damage predictions obtained with the radial mesh just beyond peak load, d = 1/8 inch, 
m = 4. 


curve (C to D), corresponding to extensive delamination between the -45° and 0° plies, prior to reaching the ultimate 
load. The predicted average failure stress is 20% higher than observed in the tests, and the predicted ultimate stress 
is approximately 10% higher than obtained in the tests (Table 4). The radial mesh, on the other hand fails to predict 
the delamination failure mechanism, and ultimate failure is predicted early. The predicted load versus displacement 
response obtained with the radial mesh is linear up to ultimate failure, and the ultimate load is approximately 25% 
lower than the ultimate load obtained in the tests. 

Matrix damage predictions in the individual plies obtained from the radial mesh at the peak load are shown in 
Fig. 10. At the peak load, some axial- splitting damage has developed in the 0° layer at the hole edge. The predicted 
matrix damage trajectories in the 45° and 90° plies are very similar, and very little delamination develops between 
these plies. Fairly significant delaminations develop at the 90/-45 interface and at the -45/0 interface (not shown). 
However, significant stress relief at the hole edge is not predicted, and the specimen is predicted to fail prematurely 
as a result of fiber failure in the 0° ply. 

Damage predictions obtained with the aligned mesh at selected points in the load-displacement response are 
provided in Fig. 11. The analysis results indicate that matrix cracking develops in all plies at low load levels (Point 
A, 3854 lb). With increasing load, significant progressive growth of delaminations begins to form beyond the hole- 
edge (Point B, 4496 lb). With further increase in load, the 45° surface cracks propagate to the free edge, and a 
delamination develops at the 45/90 interface between the laminate free edges and the 45° cracks. The edge 
delaminations then link-up with the delaminations around the edge of the hole, forming a delamination over the 
entire width of the specimen (Point C, 4802 lb). The damage is then able to propagate through the laminate 
thickness, and a major delamination between the -45° ply and the 0° ply develops (Point D, 4474 lb). The 
delamination propagates rapidly throughout the specimen, leading to nearly complete separation of the -45° and 0° 
layers and the load drop shown in Fig. 9. In addition, extensive delamination between the 90° and 45° layers, 
bounded by cracking in the off-axis plies is observed. With continued loading after the load drop, the axial splits 
propagate to the specimen ends and the interface between the -45° and 0° plies completely delaminates. At the 
ultimate stress, delamination failure of the specimen is predicted, with fiber failure in the 0-degree ply. 

Comparison of damage predictions with typical X-Ray/CT (Computed Tomography) and ultrasonic images 
taken from interrupted tests of specimens conducted by the authors are provided in Figs. 12 and 13 for a specimen 
with a 1/8-inch diameter hole, and m = 4. The test specimens were incorrectly fabricated and were unsymmetrical, 
with the -45° ply on one half of the laminate thickness thicker than on the other half, so only a qualitative 
comparison can be made. 


11 

American Institute of Aeronautics and Astronautics 






A (3854 lb) 


B (4496 lb) 


C (4802 lb) 


D (4474 lb) 








Damage 

Index 


Fully 

Damaged 



No 

Damage 


Figure 11. Damage evolution predictions obtained with the aligned mesh, d = 1/8 inch, m = 4. 
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Loading direction 



(a) X-Ray/CT, matrix cracks (b) X-Ray/CT, delaminations at 45/90 interface 



(c) Matrix crack predictions, all plies (d) Predicted delaminations at 45/90 interface 

Figure 12. Experimental damage and damage predictions near peak experimental load and peak 
numerical load, d = 1/8 inch, m = 4. 

In Fig. 12, predictions obtained near the predicted failure load are compared with NDE images of a test specimen 
interrupted at 95% of the experimentally observed failure load. The results shown in Fig. 12 indicate that the 
analysis predicts well the delamination pattern at the 45/90 interface. Additionally, the predicted matrix crack 
pattern captures the dominant cracks in the 45° surface ply, and the dominant axial splits. The axial splits observed 
in the tests develop asymmetrically around the hole, where the asymmetry in the splits is not accurately captured in 
the simulations. The experimental results show distributed matrix cracks in the -45° ply developing along the splits 
in the 45° and 0° plies as well as distributed matrix cracks in the 90° ply along the splits in the 45° ply, where the 
analysis predicts only a few cracks in these plies. In Fig. 13, X-Ray/CT images that were obtained just after the peak 
load are compared with the analysis predictions at a similar point in the load history. The analysis predictions agree 
well with the experimental delamination shapes at the 45/90 and 90/-45 interface. However, the delamination at the - 
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(a) 45/90 interface 


(b) 90/-45 interface 


(c) -45/0 interface 




(d) 45/90 interface (e) 90/-45 interface (f) -45/0 interface 

Figure 13. X-Ray/CT images of delaminations at all interfaces after first load drop and predicted 
delaminations just beyond numerical first load drop, d = 1/8 inch, m = 4. 


45/0 interface is predicted to develop both inside and outside the axial splits, while in the tests the delamination 
develops primarily outside the axial splits. 

B. In-plane and Through-Thickness Scaling Effect 

Analysis predictions obtained with the aligned mesh for studying the effect of in-plane scaling (2D scaling) and 
through-thickness scaling (ID scaling) for ply-level thickness scaled specimens are summarized and compared with 
experimental data in Table 5. For the specimen configurations considered, delamination and pull-out failure modes 
were observed in the tests. In the data presented, a load drop of 5% on the load-displacement curve is taken to 
represent specimen failure. Stress values are presented and are obtained by dividing the failure loads by the 
unnotched cross-sectional area. The results provided in the table indicate that the analysis generally predicts the 
proper failure mode. Additionally, the predicted sequence of damage development and interaction of subcritical 
damage and fiber failure agrees qualitatively with the experimental observations. Damage evolution results shown 
previously for the specimen configuration with a 1 /8-inch diameter hole, and m = 4 are representative of the 
predictions of the delamination failure mode for all such failure cases, consisting of significant matrix damage, 
delaminations between all plies, extensive delamination at the -45/0 interface, and no fiber failure. For specimens 
that were predicted to fail in the pull-out failure mode, damage development is characterized by significant sub- 
critical matrix damage in all plies, delamination at all ply interfaces, and fiber failure in the 0° ply. For equal in- 
plane dimensions, and reduced values of ply thickness (reduced value of m ), the extent of delamination at the -45/0 
interface is reduced and fiber failure in the 0° ply is more extensive. This is illustrated in Fig. 14 which compares 
damage predictions in the -45/0 interface and in the 0° ply for laminates with a 1 /8-inch diameter hole, that failed in 
the pull-out failure mode, m = 1, and in the delamination failure mode, m = 4. 
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Table 5. Test [20], Analysis Failure Stress (ksi) and Failure Mode Comparison for Open-Hole- 
Tension o f (45 fW /90 fW /-45 fW /0 fyi ) s Laminates 


Laminate Hole 

Thickness, m Diameter (in) 178 174 172 


1 

Test 

Analysis 

82.7 

64.1 

50 

2 

Test 

57.4 

72.2 


Analysis 

50.6 

49.1 

4 

Test 

39.9 

41.3 52.5 


Analysis 

49.5 

56.1 67.5 


Pull-out Delamination 

Comparison of the experimental results and the numerical predictions studying the effect of scaling on the failure 
load for specimen configurations with m = 4, and for specimen configurations with d = 1 /8-inch is also provided in 
Figs. 15a and 15b, respectively. The analyses conducted with the aligned mesh consistently predict the delamination 
failure mode for specimen configurations having m = 4. For these specimens the so-called inverse hole-size effect 
observed in the tests is predicted, as shown in Fig. 15a. The failure loads are, however, consistently over predicted. 
Experimental results for the specimens with a 1 /8-inch diameter hole and scaled in the thickness direction indicate a 
decrease in failure load with increase in ply thickness, as shown in Fig. 15b. Additionally the specimen with m = 1 
failed in the pullout failure mode, and the specimens with m = 2, 4 failed in the delamination failure mode. The 
analyses predicted failure in the pull-out failure mode for specimens with m = 1,2, and in the delamination failure 
mode for the specimen with m = 4. The failure loads for specimens exhibiting the pull-out failure are consistently 
under-predicted. 



-45/0 interface, m = 1 



-45/0 interface, m =4 
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Matrix failure, 0 C 

II 
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Matrix failure, 0° ply, m = 4 



Fiber failure, 0° ply, m = 1 


O 


Fiber failure, 0° ply, m = 4 


Figure 14. Damage predictions obtained with the aligned mesh at a point just beyond peak load, d = 
1/8 inch, m = 1 and 4. 


15 

American Institute of Aeronautics and Astronautics 




90 
80 
70 
60 

Stress, 50 
ksi 40 

30 
20 
10 
0 

Diameter, inch 



Stress, 

ksi 



Analysis 
Test [20] 


1 2 4 


Scale factor, m 


(a) In-plane scaling, for thickness m = 4 (b) Through-thickness scaling, d = 1/8 inch 

Figure 15. Experimentally observed and predicted effect of in-plane scaling and through-thickness 
scaling on failure stress. 


The critical event in determining the failure stress for the delamination failure mode is the point at which the 
surface matrix crack propagates to the laminate free edge, and the delaminations in the 45/90 interface at the hole 
edge and at the free edge link-up. In all of the analysis predictions, this event is delayed compared to the test results. 
One possible cause for the over-prediction of the failure stress is insufficient mesh refinement. Additional study is 
required to determine the effect of mesh refinement on the analysis predictions. Another possible explanation for the 
over-prediction can be obtained by comparing the experimentally determined crack patterns, just prior to the load 
drop, and the predicted crack patterns. The X-Ray/CT data provided in Fig. 12, shows multiple distributed cracks in 
the 90° ply along the split in the surface 45° ply, where the analysis predicts development of a single crack in the 
90° ply. The multiple cracks observed experimentally in the 90° ply may facilitate development of the delaminations 
at the 45/90 interface, and result in link-up of damage at the hole edge and at the free edge at a load lower than 
predicted. Additionally, thermal residual stresses from the curing process were not included in the analysis and may 
result in predicted late onset of matrix cracking and the accompanying delamination. A delayed onset and 
propagation of the sub-critical matrix and delamination damage in the analysis may also explain the under- 
prediction of the failure loads for the specimens that failed in the pull-out failure mode. In this case, the sub-critical 
damage reduces the stress concentration at the hole edge and therefore delays the onset of fiber failure observed in 
test. 


V. Concluding Remarks 

The present paper presents a study on the performance of a state-of-the-art continuum damage mechanics (CDM) 
model for intralaminar damage, coupled with cohesive elements for interlaminar damage for failure simulation of 
quasi-isotropic open-hole tension specimens. Limitations of the CDM models are discussed with an emphasis on the 
influence of element orientation on the analysis predictions. Results of these studies indicate that reliable prediction 
of matrix crack paths and stress relaxation after cracking requires a mesh with element edges aligned with the crack 
direction. Based on the findings of this study, a modeling approach is proposed for the analysis of the open-hole 
tension specimens that consists of different meshes within the individual plies, such that the element edges are 
aligned with the ply fiber direction, and cohesive element layers embedded between plies. Then, tie constraints are 
used to connect the individual ply meshes with the cohesive element layers. 

A typical radial mesh, commonly used for the analysis of open-hole specimens, and the proposed aligned mesh 
are used in an explicit analysis to simulate progressive intralaminar and interlaminar damage, and ultimate failure of 
the quasi-isotropic open-hole tension specimens. Numerical predictions are compared with experimental data 20,21 to 
assess the modeling approach for predicting sub-critical damage development, and the ultimate failure mode and 
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failure load. In the tests, three failure modes were observed for different specimen configurations: brittle, pull-out 
and delamination. The pull-out and delamination failure modes are considered in the present paper. The pull-out 
failure mode is characterized by matrix cracking in all off-axis plies, limited delamination between plies, and fiber- 
failure in the 0° plies. The delamination failure mode is characterized by matrix cracking in all off-axis plies and 
extensive delamination between plies. Comparison of the numerical predictions obtained using the radial mesh with 
the experimental data for a specimen configuration exhibiting the delamination failure mode indicates that the radial 
mesh fails to predict the delamination failure mode and the ultimate failure is predicted to occur early. The 
delamination failure mode is, however, successfully predicted with the proposed aligned mesh, and the predicted 
damage evolution qualitatively agrees with the experimentally observed damage evolution. The aligned mesh model 
is used further to simulate the affect of in-plane (2D) and through-thickness (ID) scaling on the failure mechanism 
and failure load. For the specimen configurations considered, the numerical models generally predict the proper 
failure mechanism and the proper trends in the failure load. The failure loads for specimens exhibiting the 
delamination failure load are, however, consistently over predicted, and the failure loads for specimens exhibiting 
the pull-out failure mode are consistently under-predicted. The discrepancy in strength prediction is attributed to 
predicted late onset and propagation of sub-critical damage. In the case of the delamination failure mode, late onset 
of damage results in delayed link-up of damage through the laminate width compared to the experimental results and 
over-prediction of the laminate strength. In the case of the pull-out failure mode, late onset and development of sub- 
critical damage in the analysis which results in limited reduction of the stress concentration at the hole edge, and 
early onset of fiber failure. 
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